Measuring degradation on schroeders

Fine scale acoustic perception in zebra finches

Author

Marcelo Araya-Salas

Published

August 2, 2023

Source code and data found at https://github.com/maRce10/acoustic-fine-features-zebra-finch

 

1 Purpose

  • Evaluate if schroeder structure is degraded during transmission
  • Evaluate differencial degradation in open and closed environments

 

2 Report overview

 

Code
source("~/Dropbox/Projects/geographic_call_variation_yellow-naped_amazon/scripts/MRM2.R")

warbleR_options(parallel = 22)

wave_col <- viridis(10)[7]

path <- "./data/raw/degradation_experiment"

3 Automatic annotation

Use the package baRulho to fin the position of schroeders in the re-recorded files

Code
master_annotations <- imp_raven(path = "./data/processed", files = "schroeder_master.wav.txt", all.data = TRUE, warbler.format = TRUE)


master_annotations$sound.id <- paste0("f0:", master_annotations$f0, "_comp:", master_annotations$label)

master_annotations$sound.id[1] <- "start_marker"
master_annotations$sound.id[nrow(master_annotations)] <- "end_marker"

names(master_annotations)

# table(master_annotations$sound.id)

markers_in_tests <-
    find_markers(X = master_annotations,  # annotations of sounds in master file
                 path = path, 
                    ) # path to supplementary files


aligned_tests <-
    align_test_files(X = master_annotations, # annotations of sounds in master file
                     Y = markers_in_tests, # position of markers in test files
                     path = path, # folder with files
                     output = "data.frame") 


check_sels(aligned_tests, path = path)

aligned_tests$start[grep("end_marker", aligned_tests$sound.id, invert = TRUE)] <- aligned_tests$start[grep("end_marker", aligned_tests$sound.id, invert = TRUE)] - 1.8
aligned_tests$end[grep("end_marker", aligned_tests$sound.id, invert = TRUE)] <- aligned_tests$end[grep("end_marker", aligned_tests$sound.id, invert = TRUE)] - 1.8


exp_raven(X = aligned_tests, file.name = "pl", path = path, sound.file.path = path)


aligned_tests <- aligned_tests[grep("marker", aligned_tests$sound.id, invert = TRUE), ]


aligned_tests_est <- selection_table(aligned_tests,extended = TRUE, confirm.extended = FALSE, path = path)


saveRDS(aligned_tests_est, file.path(path, "extended_sel_table_degradation_exp.RDS"))

4 Measuring schroeder dissimilarity

4.1 Dynamic-time warping pairwise distance

  • Both Schroeders have the same length
  • One is duplicated and the other one is slide across the duplicated one
  • The minimum DTW distance is kept as a dissimilarity measure
Code
mean_segment <- function(wave, cores = 1, plot = TRUE, pb = TRUE,
    thinning = 1, col = wave_col, mean = TRUE, type = "ac", npeak = 20) {
    # thin
    if (thinning < 1) {
        if (length(wave@left) * thinning < 10) {
            stop2("thinning is too high, no enough samples left for at least 1 sound file")
        }

        # reduce size of envelope
        wavefrm <- stats::approx(x = seq(0, duration(wave), length.out = length(wave@left)),
            y = wave@left, n = round(length(wave@left) * thinning),
            method = "linear")$y
    } else {
        wavefrm <- wave@left
    }

    # get empirical mode decomposition
    if (type == "EMD") {
        emds <- EMD::emd(wavefrm, seq_len(length(wavefrm)), boundary = "wave")

        perd <- emds$imf[, 4]/max(emds$imf[, 4])
        # plot(x = seq_len(length(wavefrm)), y = perd, type =
        # 'l') lines(y = wavefrm / max(wavefrm), x =
        # seq_len(length(wavefrm)), col = 'gray', lty = 2)
    }

    if (type == "ac") {
        ac <- acf(x = wavefrm, lag.max = length(wavefrm), type = "covariance",
            demean = FALSE, plot = FALSE)
        perd <- ac$acf/max(ac$acf)
    }

    tpks <- seewave::fpeaks(cbind(seq_len(length(perd)), perd), plot = FALSE,
        threshold = 0.5)

    if (nrow(tpks) > npeak) {
        tpks <- tpks[1:npeak, ]
    }

    segment_df <- data.frame(selec = seq_len(nrow(tpks)), pos = tpks[,
        1], peak = tpks[, 2])

    # get mean number of sample between peaks
    mean_dist_peak <- round(mean(diff(segment_df$pos)))

    segment_df$start <- segment_df$pos - mean_dist_peak/2
    segment_df$end <- segment_df$pos + mean_dist_peak/2

    # fix if values are out of wavefrm size
    if (segment_df$start[1] > 0) {
        segment_df$start[1] <- 0
    }
    if (segment_df$end[nrow(segment_df)] > length(wavefrm)) {
        segment_df$end[nrow(segment_df)] <- length(wavefrm)
    }

    # extract segments into a list
    segments <- lapply(seq_len(nrow(segment_df)), function(x) {
        wavefrm[segment_df$start[x]:segment_df$end[x]]
    })

    # make all the same number of samples
    segments <- lapply(segments, function(x) {
        approx(x, n = max(sapply(segments, length)))$y
    })

    # normalize between 1, -1
    segments <- lapply(segments, function(x) {
        x/max(x)
    })

    # put all segments in a data frame
    segments <- as.data.frame(segments, col.names = seq_len(length(segments)))

    # compute mean segment
    mean_segment <- rowMeans(segments)

    if (plot) {
        mean_segment_df <- data.frame(time = seq(0, 1, length.out = nrow(segments)),
            mean.amp = rowMeans(segments), sd.amp = apply(segments,
                1, sd))

        gg <- ggplot(data = mean_segment_df, mapping = aes(x = time,
            y = mean.amp)) + geom_line(color = wave_col) + geom_ribbon(aes(ymin = mean.amp -
            sd.amp, ymax = mean.amp + sd.amp), alpha = 0.2) + theme_classic(base_size = 25)

        print(gg)
    }
    if (mean) {
        return(mean_segment)
    } else {
        return(segments)
    }
}
Code
est_schr <- readRDS(file.path(path, "extended_sel_table_degradation_exp.RDS"))

mean_schroeders <- warbleR:::pblapply_wrblr_int(cl = 20, seq_len(nrow(est_schr)),
    function(x) {
        wave <- read_wave(est_schr, index = x)
        seg <- try_na(mean_segment(wave, plot = FALSE, mean = FALSE,
            type = "ac", thinning = 0.8))

        return(seg)
    })

names(mean_schroeders) <- paste0(ifelse(grepl("inside", est_schr$sound.files),
    "in_", "out_"), est_schr$sound.id)

saveRDS(mean_schroeders, file.path(path, "mean_schroeders_degradation_experiment.RDS"))
Code
mean_schroeders <- readRDS(file.path(path, "mean_schroeders_degradation_experiment.RDS"))

mean_schroeders <- mean_schroeders[!sapply(mean_schroeders, function(x) is.na(x[[1]][1]))]

mean_schroeders_list <- lapply(seq_len(length(mean_schroeders)), function(x) {
    data.frame(schroeder = names(mean_schroeders)[x], time = seq(0,
        1, length.out = nrow(mean_schroeders[[x]])), mean.amp = rowMeans(mean_schroeders[[x]]),
        sd.amp = apply(mean_schroeders[[x]], 1, sd))
})

mean_schroeders_df <- do.call(rbind, mean_schroeders_list)

4.2 Mean schroeders inside playback

  • mean period +/- standard deviation using autocorrelation
  • first 40 schroeders are shown
Code
mean_schroeders_in_df <- mean_schroeders_df[grep("in_", mean_schroeders_df$schroeder),
    ]

ggplot(data = mean_schroeders_in_df[mean_schroeders_in_df$schroeder %in%
    unique(mean_schroeders_in_df$schroeder)[1:40], ], mapping = aes(x = time,
    y = mean.amp)) + geom_line(color = wave_col) + geom_ribbon(aes(ymin = mean.amp -
    sd.amp, ymax = mean.amp + sd.amp), alpha = 0.2) + theme_classic(base_size = 5) +
    facet_wrap("~ schroeder", ncol = 5, scales = "free_y")

Code
mean_schroeders_in <- mean_schroeders[grep("in_", names(mean_schroeders))]

nms <- names(mean_schroeders_in)

cmbs <- t(combn(nms, 2))


min_dist_l <- pbapply::pbsapply(cl = 22, 1:nrow(cmbs), function(x) {
    s1 <- rowMeans(mean_schroeders_in[[cmbs[x, 1]]])
    s2 <- rowMeans(mean_schroeders_in[[cmbs[x, 2]]])

    # make same length if (length(s1) != length(s2))
    s1 <- approx(s1, n = 100)$y
    s2 <- approx(s2, n = 100)$y

    # duplicate 1
    s1 <- rep(s1, 2)

    # run dtw over longer vector
    dists <- vapply(seq_len(length(s1) - length(s2)), function(x) {
        segment <- s1[x:min(c(x + length(s2) - 1), length(s1))]

        dtw_dist <- warbleR::try_na(dtw::dtwDist(mx = rbind(s2, segment)))
        return(dtw_dist[1, 2])
    }, FUN.VALUE = numeric(1))

    return(data.frame(schr1 = cmbs[x, 1], schdr2 = cmbs[x, 2], min(dists)))
})

min_dists <- do.call(rbind, min_dist_l)

min_dists <- as.data.frame(matrix(min_dists[, 1], ncol = 3, byrow = TRUE))

names(min_dists) <- c("schr1", "schr2", "dist")

min_dists$dist <- as.numeric(min_dists$dist)

saveRDS(min_dists, file.path(path, "dtw_distance_schroeders_degradation_experiment_inside.RDS"))

4.3 Mean schroeders outside playback

Code
mean_schroeders_out <- mean_schroeders_df[grep("out_", mean_schroeders_df$schroeder),
    ]

ggplot(data = mean_schroeders_out[mean_schroeders_out$schroeder %in%
    unique(mean_schroeders_out$schroeder)[1:40], ], mapping = aes(x = time,
    y = mean.amp)) + geom_line(color = wave_col) + geom_ribbon(aes(ymin = mean.amp -
    sd.amp, ymax = mean.amp + sd.amp), alpha = 0.2) + theme_classic(base_size = 5) +
    facet_wrap("~ schroeder", ncol = 5, scales = "free_y")

Code
mean_schroeders_out <- mean_schroeders[grep("out_", names(mean_schroeders))]

nms <- names(mean_schroeders_out)

cmbs <- t(combn(nms, 2))


min_dist_l <- pbapply::pbsapply(cl = 22, 1:nrow(cmbs), function(x) {
    s1 <- rowMeans(mean_schroeders_out[[cmbs[x, 1]]])
    s2 <- rowMeans(mean_schroeders_out[[cmbs[x, 2]]])

    # make same length if (length(s1) != length(s2))
    s1 <- approx(s1, n = 100)$y
    s2 <- approx(s2, n = 100)$y

    # duplicate 1
    s1 <- rep(s1, 2)

    # run dtw over longer vector
    dists <- vapply(seq_len(length(s1) - length(s2)), function(x) {
        segment <- s1[x:min(c(x + length(s2) - 1), length(s1))]

        dtw_dist <- warbleR::try_na(dtw::dtwDist(mx = rbind(s2, segment)))
        return(dtw_dist[1, 2])
    }, FUN.VALUE = numeric(1))

    return(data.frame(schr1 = cmbs[x, 1], schdr2 = cmbs[x, 2], min(dists)))
})

min_dists <- do.call(rbind, min_dist_l)

min_dists <- as.data.frame(matrix(min_dists[, 1], ncol = 3, byrow = TRUE))

names(min_dists) <- c("schr1", "schr2", "dist")

min_dists$dist <- as.numeric(min_dists$dist)

saveRDS(min_dists, file.path(path, "dtw_distance_schroeders_degradation_experiment_outside.RDS"))

5 Statistical analysis

Modeling: - Multiple Regression on distance Matrices - Model:
\[\begin{align*} Dissimilarity &\sim frequency + components + sign \end{align*}\] - Response values scaled to make effect sizes comparable across models - Predictors were coded as pairwise binary matrices in which 0 means that calls in a dyad belong to the same level and 1 means calls belong to different levels - One model for each environment treatment as well as on the original model sounds

5.1 Original (model) Schroeders

Code
(dtw_wv_mod <- readRDS("./data/processed/matrix_correlation_dtw_distance_random_start.RDS"))
$coef
             dtw_dist  pval
Int        -1.5156699 1e-04
frequency   0.1859515 1e-04
components  1.0638637 1e-04
sign        0.6755587 1e-04

$r.squared
       R2      pval 
0.1624519 0.0001000 

$F.test
        F    F.pval 
2803.4517    0.0001 

$n
[1] 295

5.2 Outside playback

Code
min_dists <- readRDS(file.path(path, "dtw_distance_schroeders_degradation_experiment_outside.RDS"))

min_dists$dist <- scale(min_dists$dist)

dist_tri <- PhenotypeSpace::rectangular_to_triangular(min_dists)

freq_bi_tri <- as.dist(binary_triangular_matrix(group = sapply(strsplit(rownames(dist_tri),
    "_"), "[[", 2)))

comp_bi_tri <- as.dist(binary_triangular_matrix(group = gsub("_n|_p",
    "", sapply(strsplit(rownames(dist_tri), "_"), "[[", 3))))

sign_bi_tri <- as.dist(binary_triangular_matrix(group = sapply(strsplit(rownames(dist_tri),
    "_"), "[[", 4)))

rect_var <- cbind(as.dist(dist_tri), freq_bi_tri, comp_bi_tri, sign_bi_tri)

colnames(rect_var) <- c("dtw_dist", "frequency", "components", "sign")

dtw_wv_mod <- MRM2(formula = dtw_dist ~ frequency + components + sign,
    nperm = 10000, mat = rect_var)

saveRDS(dtw_wv_mod, "./data/processed/matrix_correlation_dtw_distance_outside_experiment.RDS")
Code
(dtw_wv_out_mod <- readRDS("./data/processed/matrix_correlation_dtw_distance_outside_experiment.RDS"))
$coef
               dtw_dist   pval
Int        -0.017062308 0.0001
frequency   0.010099084 0.0001
components  0.001829854 0.6014
sign        0.013210074 0.0001

$r.squared
          R2         pval 
5.547013e-05 1.000000e-04 

$F.test
        F    F.pval 
0.7641249 0.0001000 

$n
[1] 288

5.3 Inside playback

Code
min_dists <- readRDS(file.path(path, "dtw_distance_schroeders_degradation_experiment_inside.RDS"))

min_dists$dist <- scale(min_dists$dist)

dist_tri <- PhenotypeSpace::rectangular_to_triangular(min_dists)

freq_bi_tri <- as.dist(binary_triangular_matrix(group = sapply(strsplit(rownames(dist_tri),
    "_"), "[[", 2)))

comp_bi_tri <- as.dist(binary_triangular_matrix(group = gsub("_n|_p",
    "", sapply(strsplit(rownames(dist_tri), "_"), "[[", 3))))

sign_bi_tri <- as.dist(binary_triangular_matrix(group = sapply(strsplit(rownames(dist_tri),
    "_"), "[[", 4)))

rect_var <- cbind(as.dist(dist_tri), freq_bi_tri, comp_bi_tri, sign_bi_tri)

colnames(rect_var) <- c("dtw_dist", "frequency", "components", "sign")

dtw_wv_mod <- MRM2(formula = dtw_dist ~ frequency + components + sign,
    nperm = 10000, mat = rect_var)

saveRDS(dtw_wv_mod, "./data/processed/matrix_correlation_dtw_distance_inside_experiment.RDS")
Code
(dtw_wv_in_mod <- readRDS("./data/processed/matrix_correlation_dtw_distance_inside_experiment.RDS"))
$coef
             dtw_dist  pval
Int        -1.0128809 1e-04
frequency   0.3789239 1e-04
components  0.2454399 1e-04
sign        0.9017830 1e-04

$r.squared
       R2      pval 
0.2208956 0.0001000 

$F.test
        F    F.pval 
4070.1899    0.0001 

$n
[1] 294

5.4 Combined results

Code
mods <- list(dtw_wv_mod = dtw_wv_mod, dtw_wv_in_mod = dtw_wv_in_mod,
    dtw_wv_out_mod = dtw_wv_out_mod)

names(mods) <- c("Original Schroeders", "Inside playback", "Outside playback")

estimates <- do.call(rbind, lapply(seq_along(mods), function(x) {
    Y <- data.frame(mods[[x]]$coef[-1, ])
    Y$rel_coef <- Y[, 1]/max(Y[, 1])
    Y$mod <- names(mods)[x]
    Y$predictor <- rownames(Y)
    names(Y) <- c("coef", "p", "rel_coef", "model", "predictor")
    return(Y)
}))

estimates$rel_coef <- ifelse(estimates$p < 0.05, estimates$rel_coef,
    0)
estimates$signif <- ifelse(estimates$p < 0.05, "p < 0.05", "p >= 0.05")

estimates$model <- factor(estimates$model, levels = c("Outside playback",
    "Inside playback", "Original Schroeders"))

ggplot(estimates, aes(x = predictor, y = model, fill = rel_coef)) +
    geom_tile() + coord_equal() + scale_fill_gradient2(low = viridis(10)[3],
    high = viridis(10)[7], guide = "none") + geom_text(aes(label = round(coef,
    3), color = signif)) + scale_color_manual(values = c("black",
    "gray")) + labs(x = "", y = "", color = "P value") + theme_classic() +
    theme(axis.text.x = element_text(color = "black", size = 11, angle = 30,
        vjust = 0.8, hjust = 0.8))

Effect sizes per model and predictor (color intensity shows effect size magnitude relative to the highest effect size within the model).

 

 

6 Takeaways

  • Schroeder structure fairly detectable at ~30 cm
  • Schroeder discrimination mostly affected in outside playback
  • Harmonic content seems to be the most affect feature, both inside and outside

 

Session information

R version 4.1.0 (2021-05-18)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.2 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/atlas/libblas.so.3.10.3
LAPACK: /usr/lib/x86_64-linux-gnu/atlas/liblapack.so.3.10.3

locale:
 [1] LC_CTYPE=pt_BR.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=es_CR.UTF-8        LC_COLLATE=pt_BR.UTF-8    
 [5] LC_MONETARY=es_CR.UTF-8    LC_MESSAGES=pt_BR.UTF-8   
 [7] LC_PAPER=es_CR.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=es_CR.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] numform_0.7.0        ecodist_2.0.9        PhenotypeSpace_0.1.0
 [4] ggplot2_3.4.2        baRulho_1.1.0        ohun_1.0.0          
 [7] Rraven_1.0.13        viridis_0.6.3        viridisLite_0.4.2   
[10] warbleR_1.1.28       NatureSounds_1.0.4   tuneR_1.4.4         
[13] seewave_2.2.1        formatR_1.11         knitr_1.43          
[16] kableExtra_1.3.4    

loaded via a namespace (and not attached):
 [1] colorspace_2.1-0      rjson_0.2.21          deldir_0.2-10        
 [4] class_7.3-22          xaringanExtra_0.7.0   rstudioapi_0.14      
 [7] proxy_0.4-27          spatstat.data_2.1-0   farver_2.1.1         
[10] Deriv_4.1.3           remotes_2.4.2         fansi_1.0.4          
[13] xml2_1.3.5            codetools_0.2-18      splines_4.1.0        
[16] polyclip_1.10-0       jsonlite_1.8.7        vdiffr_1.0.5         
[19] packrat_0.9.0         cluster_2.1.2         png_0.1-7            
[22] spatstat.sparse_2.0-0 compiler_4.1.0        httr_1.4.4           
[25] backports_1.4.1       assertthat_0.2.1      Matrix_1.5-4.1       
[28] fastmap_1.1.1         cli_3.6.1             htmltools_0.5.5      
[31] tools_4.1.0           igraph_1.5.0.1        gtable_0.3.3         
[34] glue_1.6.2            dplyr_1.0.10          Rcpp_1.0.11          
[37] raster_3.4-13         vctrs_0.6.3           svglite_2.1.0        
[40] nlme_3.1-162          xfun_0.39             stringr_1.5.0        
[43] brio_1.1.3            testthat_3.1.10       rvest_1.0.3          
[46] lifecycle_1.0.3       goftest_1.2-2         MASS_7.3-60          
[49] scales_1.2.1          spatstat.core_2.3-0   spatstat.utils_2.2-0 
[52] parallel_4.1.0        Sim.DiffProc_4.8      yaml_2.3.7           
[55] pbapply_1.7-2         gridExtra_2.3         rpart_4.1-15         
[58] stringi_1.7.12        e1071_1.7-13          checkmate_2.2.0      
[61] permute_0.9-5         rlang_1.1.1           pkgconfig_2.0.3      
[64] systemfonts_1.0.4     dtw_1.23-1            bitops_1.0-7         
[67] evaluate_0.21         lattice_0.21-8        tensor_1.5           
[70] sf_1.0-14             labeling_0.4.2        htmlwidgets_1.5.4    
[73] tidyselect_1.2.0      magrittr_2.0.3        R6_2.5.1             
[76] fftw_1.0-7            generics_0.1.3        DBI_1.1.3            
[79] pillar_1.9.0          withr_2.5.0           mgcv_1.8-42          
[82] units_0.8-2           abind_1.4-5           RCurl_1.98-1.12      
[85] sp_1.5-1              tibble_3.2.1          crayon_1.5.2         
[88] KernSmooth_2.23-21    utf8_1.2.3            spatstat.geom_2.2-2  
[91] rmarkdown_2.23        grid_4.1.0            sketchy_1.0.2        
[94] vegan_2.5-7           digest_0.6.33         classInt_0.4-9       
[97] webshot_0.5.4         signal_0.7-7          munsell_0.5.0